本篇博客的目标：对数似然概率公式的Matlab编程

对数似然概率公式：

$\begin{aligned} \ln p(\mathbf{X} | \boldsymbol{\mu}, \mathbf{\Sigma})= -\frac{1}{2}(Dln(2\pi)+ln|\Sigma|+\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)^{T} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)) \end{aligned}$

编程难点：公式中的$x$是列向量，而Matlab中的$X$是矩阵，为了提高代码的效率避免使用过多的$for$语句，这里使用向量化技术完成算法的编程任务。

公式推导

参考链接：https://www.cnblogs.com/ccienfall/p/6049021.html

多维高斯函数的似然概率：

$\mathcal{N}(\mathbf{x} | \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T} \mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}$

若是考虑数据集$\mathbf{X}=\left(\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}\right)^{T}$，再求对数则为：

$\ln p(\mathbf{X} | \boldsymbol{\mu}, \mathbf{\Sigma})=-\frac{N D}{2} \ln (2 \pi)-\frac{N}{2} \ln |\mathbf{\Sigma}|-\frac{1}{2} \sum_{n=1}^{N}\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)^{T} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)$

⭐️但是我需要的似然概率，$x_i$在第$k$个分布下的似然概率：【实数】

$\mathcal{N}(\mathbf{x_i} | \boldsymbol{\mu_k}, \mathbf{\Sigma_k})=\frac{1}{(2 \pi)^{D / 2}} \frac{1}{|\mathbf{\Sigma_k}|^{1 / 2}} \exp \left\{-\frac{1}{2}(\mathbf{x_i}-\boldsymbol{\mu_k})^{T} \mathbf{\Sigma_k}^{-1}(\mathbf{x_i}-\boldsymbol{\mu_k})\right\}$

展开后，得证：

$\begin{aligned} \ln p(\mathbf{X} | \boldsymbol{\mu}, \mathbf{\Sigma})= &-\frac{ D}{2} \ln (2 \pi)-\frac{1}{2} \ln |\mathbf{\Sigma}|-\frac{1}{2} \left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)^{T} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right) \\ = & -\frac{1}{2}(Dln(2\pi)+ln|\Sigma|+\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)^{T} \mathbf{\Sigma}^{-1}\left(\mathbf{x}_{\mathbf{n}}-\boldsymbol{\mu}\right)) \end{aligned}$

最终代码

function logprob = lgmmprob(data, mu, sigma, w)
% 计算GMM的观察概率（似然概率）
ndim = size(data, 1);
C = sum(mu.*mu./sigma) + sum(log(sigma));
D = (1./sigma)' * (data .* data) - 2 * (mu./sigma)' * data  + ndim * log(2 * pi);
logprob = -0.5 * (bsxfun(@plus, C',  D));
logprob = bsxfun(@plus, logprob, log(w));

舍弃系数-0.5，进一步简化需要编程的部分：

$Dln(2\pi)+ln|\Sigma|+\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right)^{T} \mathbf{\Sigma_k}^{-1}\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right)$

需要说明的是上面所有的小项都是实数。

向量化技术⭐️

1. 改造公式

$\Sigma$在数学推导中的公式中是 $full -variance$ 矩阵，但是一般代码中的协方差矩阵默认是对角协方差矩阵。故这里就存在一个算法和实际编程中的一个不匹配的地方，故这个时候需要对原本的算法公式进行改造

• 正常的推到方程

$\left(a_{1}, a_{2}\right)*\left(\begin{array}{cc} {\Delta_{1}} & {0} \\ {0} & {\Delta_{2}} \end{array}\right)*\left(\begin{array}{l} {a_{1}} \\ {a_{2}} \end{array}\right)=\Delta_{1} a_{1}^{2}+\Delta_{2} a_{2}^{2}$

• 代码中待编程的方程

$\left(\Delta_{1}, \Delta_{2}\right)*\left(\begin{array}{l} {a_{1}} \\ {a_{2}} \end{array}\right).*\left(\begin{array}{l} {a_{1}} \\ {a_{2}} \end{array}\right)=\Delta_{1} a_{1}^{2}+\Delta_{2} a_{2}^{2}$

故原先的第三项为：$\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right)^{T} \mathbf{\Sigma_k}^{-1}\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right)$

其中令$\Sigma$是一个列向量，可以见下节的变量维度说明。$(\Delta_1,\Delta_2)=(\Sigma_{k}^{-1})^T$

$(\mathbf{\Sigma_k}^{-1})^T*\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right).*\left(\mathbf{x}_i-\boldsymbol{\mu_k}\right)$

• 本例中，代码所需编写的方程为：

$Dln(2\pi)+ln|\Sigma_k|+(\mathbf{\Sigma_k}^{-1})^T*[x_i\ .*\ x_i-2\mu_k\ .*\ x_i+\mu_k\ .*\ \mu_k]$

2. 变量维度说明：

$x_i:(p,1)$

$X=\{x_1,x_2,\dots,x_i,\dots,x_n\}:(p,n)$

$\Sigma_k:(p,1)$

$\Sigma = \{\Sigma_1,\Sigma_2,\dots,\Sigma_k,\dots,\Sigma_K\}:(p,K)$

$\mu_k:(p,1)$

$\mu=\{\mu_1,\mu_2,\dots,\mu_k,\dots,\mu_K\}:(p,K)$

3. 列向量改矩阵⭐️

这个是最关键的一步，一步一步解答，以便下次遇到的时候可以轻松应对，上式是原数学方程式中的数学向量，但为了编程，需要向量化。

• 首先明确向量化的目的：用$\Sigma,\mu,X$代替$\Sigma_k,\mu_k,x_i$
• 其次明确维度：每一项的都是实数

这里再明确一下待编程的方程为：

$Dln(2\pi)+ln|\Sigma_k|+(\mathbf{\Sigma_k}^{-1})^T*[x_i\ .*\ x_i-2\mu_k\ .*\ x_i+\mu_k\ .*\ \mu_k]$

根据每一小项是否与分布$k$和样本$i$有关可进一步划分：

	与k有关	与k无关
与$i$有关	1. k与i 矩阵乘 2. k 与 i 点乘，需转换为矩阵乘 3. 不满足以上两个条件
与$i$无关	$ln\|\Sigma_k^{-1}\|$	$Dln(2\pi)$

1. 情况一：与$k$和$i$都无关的项，$Dln(2\pi):$实数

ndim * log(2 * pi)

2. 情况二：仅与$k$有关的项，$ln|\Sigma_k^{-1}|:$ 实数

   这里试求方差的行列式，又方差是对角矩阵，所以可简化为：

$ln|\left(\begin{array}{cc} {\Delta_{1}} & {0} \\ {0} & {\Delta_{2}} \end{array}\right)|=ln|\Delta_1\cdot\Delta_2|=ln\Delta_1+ln\Delta_2=sum(ln(\Delta_1,\Delta_2))$

​ $ln|\Sigma_k|=sum(\Sigma_k):(1,1)$：

​ 可以向右扩展：

​ $\{sum(\Sigma_1),sum(\Sigma_2),\dots,sum(\Sigma_k)\}=sum(\Sigma):(1,k)$

3. 情况三：$k$与$i$已经分开，且由矩阵乘相连

   $(\mathbf{\Sigma_k}^{-1})^T*[x_i\ .*\ x_i]$：$(p,1)^T*(p,1).*(p,1)=(1,1)$实数

   向量化具体细节讲解：

   元素向量化：将$x_i$横向扩展，将$\Sigma_k$纵向扩展

$\{(\mathbf{\Sigma_k}^{-1})^Tx_1.*x_1,(\mathbf{\Sigma_k}^{-1})^Tx_2.*x_2,\dots,(\mathbf{\Sigma_k}^{-1})^Tx_n.*x_n\} \\ = (\mathbf{\Sigma_k}^{-1})^T *\{x_1.*x_1,\ x_2.*x_2,\ \dots,x_n.*x_n\} \\ =(\mathbf{\Sigma_k}^{-1})^T*(X.*X):(1,p)*(p,n)=(1,n)还可以向上下扩展\\ =((\mathbf{\Sigma_1}^{-1},\mathbf{\Sigma_2}^{-1},\dots,\mathbf{\Sigma_K}^{-1}))^T*(X.*X) \\=(\Sigma^{-1})^T*(X.*X):(k,n)$

4. 情况四：$k$和$i$为点乘，需要转换为矩阵乘：$(\mathbf{\Sigma_k}^{-1})^T*[-2\mu_k.*x_i]$：$(p,1)^T*(p,1)=(1,1)$实数

但此时有个问题，无法进行向量化，因为这里的运算规则必须先进行点乘运算，再进行矩阵乘运算，若需向量化中间必须用矩阵乘相连，原因后面解释。

$-2*(\Sigma_k^{-1})^T.*\mu_k^T*x_i=-2(\Sigma_k^{-1}.*\mu_k)^T*x_i$

​ 向量化处理同情况三：直接去掉下标即可

$-2(\Sigma^{-1}.*\mu)^T*X$

5. 情况五：都为$k$，中间既有矩阵乘又有点乘

   $(\mathbf{\Sigma_k}^{-1})^T*\mu_k\ .*\ \mu_k:(p,1)^T(p,1)=(1,1)$

   将矩阵乘转换为点乘：

   $(\mathbf{\Sigma_k}^{-1})^T*\mu_k\ .*\ \mu_k=sum(\Sigma_k^{-1}.*\mu_k.*\mu_k)$

对于一些细节的补充：

1. 为啥必须为矩阵乘来连接两个部分,举例如下：

   若：$X=(x_1,x_2,\dots,x_n)$

   $(Ax_1,Ax_2,\dots,Ax_n)=A*(x_1,x_2,\dots,x_n)=A*X$

   而点乘则不满足这条性质：

   $(A.*x_1,A.*x_2,\dots,A.*x_n)*(x_1,x_2,\dots,x_n)=A*X$

2. 为啥情况五中，不允许出现点乘？

   $fifth=\{(\mathbf{\Sigma_1}^{-1})^T*(\mu_.*\mu_1),(\mathbf{\Sigma_2}^{-1})^T*(\mu_2.*\mu_2),\dots,(\mathbf{\Sigma_k}^{-1})^T*(\mu_k.*\mu_k)\}$

   根据向量化的定理，根本没办法凑出一个完整的矩阵如$\mu$，若直接去除下标则为：

   $(\Sigma^{-1})^T*\mu.*\mu:(p,k)^T(p,k)=(k,k)$

   上式肯定不可能办到，因为向量化同一个小标只能往一个方向扩展，即应为(1,k)